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Abstract 

We investigate the effect of electric synapses (gap junctions) on col- 
lective neuronal dynamics and spike statistics in a conductance-based 
Integrate-and-Fire neural network, driven by a Brownian noise, where 
conductances depend upon spike history. We compute explicitly the time 
evolution operator and show that, given the spike-history of the network 
and the membrane potentials at a given time, the further dynamical evo- 
lution can be written in a closed form. We show that spike train statistics 
is described by a Gibbs distribution whose potential can be approximated 
with an explicit formula, when the noise is weak. This potential form 
encompasses existing models for spike trains statistics analysis such as 
maximum entropy models or Generalized Linear Models (GLM). We also 
discuss the different types of correlations: those induced by a shared stim- 
ulus and those induced by neurons interactions. 

Keywords Neural networks dynamics; spike train statistics; Gibbs distri- 
butions. 

1 Introduction 

Communication between neurons involves chemical and electric synapses. Elec- 
tric synapses transmission is mediated by gap junctions with a direct electrical 
communication between cells ^21) , allowing faster communication than chemical 
synapses. Electrical coupling between cells can be found in many parts of the 
nervous system [3], [11], and also outside: for example, certain cells in the heart 
and pancreas are connected by gap junctions [Hj. At the network level, elec- 
tric synapses have several prominent effects such as neurons synchronization [2] , 
[26j . |63j . and the generation of neural rhythms [37], [3]. In the retina, gap junc- 
tions are involved in circuits containing the five neuronal types (photo-receptors, 
horizontal, bipolar, amacrine and ganglion cells). In particular, amacrine and 
ganglion cells express numerous gap junctions, resulting in extensive electrical 
networks of low-resistance pathways for direct electrical signaling between cells 
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[55] , Unlike studies carried out in other parts of the Central Nervous System 
that often assign generic functions to electrical coupling such as increased spike 
synchrony, studies in the retina have been able to detail the specific roles of 
individual gap junctions in visual processing. Several distinct functional roles 
have been elucidated for retinal gap junctions, related not only to the propa- 
gation of signals, but also to the encoding of specific visual information [3]. In 
particular, electric synapses play a role in detection of approaching motion, and 
detection of differential motion f37 . 

Dealing with spike population coding, it is clear that neurons interactions 
highly constrain the collective spike response of a neural assembly to stimuli. 
Therefore, unveiling the respective effect of chemical and electric synapses on 
spike responses is a mandatory step toward a better understanding of spike 
coding. As an example, the retina encodes visual stimuli in a highly parallel, 
adaptive and computationally efficient way, extracting specific stimulus features 
in a highly selective manner. Central to this ability is the segregation of reti- 
nal responses into a large number of parallel pathways, resulting in a mosaic of 
functional units encoding different aspects of the visual stimulus emerging from 
dedicated local retinal circuits (for a review see [31]). Those circuits involve 
both electric and chemical synapses. 

On theoretical grounds, the role of gap junctions in shaping collective dynam- 
ics has been quite less studied than the role of chemical synapses, although differ- 
ent models and approaches have been used to address this problem, in the con- 
text of pattern formation, using techniques such as: Poincare map PTI ISUl [S5] . 
Lyapunov functions j57j , mean-field approach and Fokker-Planck equation [63] , 
variance analysis [56], and phase plane analysis [20l [21]. The effects of gap 
junctions on spike trains statistics is even less known. 

The goal of this paper is to push one step forward the mathematical analysis 
of the join effects of chemical and electric synapses on neurons dynamics and 
spike statistics. We propose and analyze a conductance-based Integrate-and- 
Fire model submitted to noise, inspired from the paper [74], with an extension 
including gap junctions. This work follows the paper [T3] where only the effects 
of chemical synapses were studied. It was especially shown that spike statis- 
tics is characterized by a Gibbs distribution whose potential can be explicitly 
computed. This provided moreover a firm theoretical ground for recent studies 
attempting to describe experimental rasters in the retina |88j as well as in the 
parietal cat cortex [5 3) by Gibbs distributions and maximum entropy principle. 
In the present paper, we extend the mathematical analysis of |13j including 
electric synapses, allowing us to consider the simultaneous effects of both types 
of synapses in a single formalism. 

The main advantage of this type of model, rendering the mathematical anal- 
ysis tractable, is that the sub-threshold dynamics of membrane potential is de- 
scribed by a linear system of stochastic differential equations (SDE). This system 
is nevertheless highly non trivial, as we show, since it is non autonomous, non 
homogeneous, with, additionally, a flow depending on the whole (spike) his- 
tory via the chemical synapses. Moreover, electric synapses introduce another 
mechanism of history dependence, where the membrane potential of a neuron 
depends on its past values, even those which are anterior to the last firing and 
reset of that neuron. The flow of this linear system of stochastic differential 
equations can be explicitly written. 
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From this, one can compute a family of transition probabilities, character- 
izing the probability of a spike pattern at a given time, given the past. These 
transition probabilities define a Gibbs distribution characterizing spike train 
statistics. The potential of this Gibbs distribution can be approached by an 
explicit form, as we show. We have therefore an explicit characterization of 
spike statistics in a model including chemical and electric synapses. 

This result has several implications in the realm of biological spike trains 
analysis. Especially, as we show, the Gibbs potential of our model encompasses 
existing models for spike trains statistics analysis such as maximum entropy 
models [7i[Sl[Hni[HIl[SIl[Ml[2SllHH] and Generalized-Linear Models (GLM) 
[711 [Sni m [Zni m] • Moreover, our formalism affords the study non stationary dy- 
namics, while stationarity is a major assumption when using maximum entropy 
models. Additionally, as we show, gap junctions introduce major correlations 
effects ruining the hope of having conditionally-upon-the-past independent neu- 
rons, a central hypothesis in GLM. Finally, three types of effects are responsible 
for neuron correlations (pairwise and higher order): shared stimulus, chemical 
couplings, and electric couplings. The two last types of correlations persists 
even when the stimulus is switched-off. 

The paper is organized as follows. In section [2] we present the model, and 
detail how to handle properly the continuous-time dynamics of membrane po- 
tentials and the discrete-time dynamics of spikes. In section [3] we give explicit 
solutions for the system of equations presented in section [2] where the network 
parameters as well as the network spike history dependence is explicit. In sec- 
tion |4] we describe the probability of spike events and the corresponding Gibbs 
distribution. The section |5] is devoted to analyzing the possible consequences 
of those mathematical results on current neuroscience studies in the domain of 
MEA (Multi Electrodes Array) retina spike train analysis. Several conclusions 
and perspectives are drawn in the last section. 

2 Model definition 

The Integrate-and-Fire model remains one of the most ubiquitous model for 
simulating and analyzing the dynamics of neuronal circuits. Despite its sim- 
plified nature, it captures some of the essential features of neuronal dynamics 
(see [49l SHI [3 [9l |47] for a review) . In this work we consider an extension of 
the conductance based Integrate-and-Fire neuron model proposed in [71]. The 
mo del- definition follows the presentation given in [TH [T3] . Dynamics is ruled 
by a set of stochastic differential equations where parameters, corresponding to 
chemical conductances, depend upon a sequence of discrete variables summa- 
rizing the action potentials emitted in the past by the neurons. In this way, the 
dynamical system defined here is ruled both by continuous time and discrete 
time dynamical variables. Let us first introduce the discrete time variables. 

2.1 Spike trains 

We consider a network of N neurons. We assume that there is a minimal time 
scale, S, set to 1 without loss of generality such that a neuron can at most emit 
an action potential (spike) within a time window of size S. This provides a time 
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discretization labeled with an integer time n. To each neuron k and discrete 
time n one associates a spike variable ujk {n) = 1 if neuron k as emitted a spike in 
the time window [nS, (n + 1)S[ and Wfe(n) = otherwise. The spike-state of the 

entire network in time bin n is thus described by a vector Lu{n) [^kin)]^_i, 
called a spiking pattern. A spike block is a finite ordered list of such vectors, 
written: 

where spike times have been prescribed between time ni to ri2- A spike train 
is a spike block i^^j^. Obviously, real spike trains start from some initial time 
to > —oo and end at some final time T < +00, but, on mathematical grounds 
the consideration of bi-infinite sequences simplifies the analysis. To alleviate 
notations we simply write lu for a spike train. 

The set of spiking pattern is denoted ^ = { 0, 1 }^ . The set of spike blocks 
is denoted A^. The set of spike trains is denoted = A^. 
To each raster uj and each neuron index j = {1 . . . N} we associate an ordered 
(generically infinite) list of "firing times" {t^^\ui)}^^ such that t''j \oj) is the 
r-th time of firing of neuron j in the raster oj. Tk{t,u)) denote the last firing 
time of neuron k. 

We note [t] the largest integer which is<t (thus [— f -2] = —2 and [1.2] = 1) 
For a function / : R x — >■ R, when we write f{t,uj) we mean /(t, wL*]^)- This 
corresponds to the physical notion of causality. The fimctions f{t,u)) that we 
consider depend, at time t, upon spikes occurring before t. 



2.2 Membrane potential dynamics 

Neurons are considered here as points, with neither spatial extension nor bio- 
physical structure (axon, soma, dendrites) . We note Vk (t) the membrane poten- 
tial of neuron fc = 1, . . . , A^, at time t. Denote V{t) the vector with entries Vi~{t). 
The continuous-time dynamics of V{t) is defined as follows. Fix a real variable 
^ > called "firing threshold" . For a fixed time t, we have two possibilities: 

1. Either Vk{t) < 6, Vfc = 1,...,A''. This corresponds to sub-threshold dy- 
namics. 

2. Or, 3k, Vk{t) > 6. Then, we speak of firing dynamics. 
2.2.1 Sub-threshold dynamics 

The sub-threshold variation of the membrane potential of neuron k at time t is 
given by: 

= -gL,k{yk-EL)-Y,gkj{t,Lo){Vk-Ej)-Y,W] (Vk -Vj)+Ik{t). (1) 
j 3 

• Ck is the membrane capacity of neuron k. 

• The term Ik{t) is a current given by: 

Ik{t) = it'\t) + aBm, (2) 
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where i^/^'^^Ht) is a deterministic external current ("stimulus"). The noise 
term ^k{t) mimics: (i) the random variation in the ionic flux of charges 
crossing the membrane per unit time at the post synaptic button, upon 
opening of ionic channels due to the binding of neurotransmitter, (ii) the 
fluctuations in the current resulting from the large number of opening and 
closing of ion channels |781 146] : (iii) noise coming from electrical synapses. 
It is common to model this noise by a Wiener white noise (diffusion ap- 
proximation). We use this modelling choice in the paper. We note <7b > 
the amplitude (mean-square deviation) of S,k{t). 

• gL,k is a leak conductance. El < is the leak reversal potential. 

• gkj{t,u)) corresponds to the conductance of the chemical synapses from 
the pre-synaptic neuron j to the post-synaptic neuron fc, while Ej is the 
reversal potential associated with the neurotransmitter emitted by neu- 
ron j. In this model, the conductance gkj{t,Lo) at time t depends upon 
the spikes emitted by the pre-synaptic neuron j. The description of this 
dependence is made explicit in section [2.2.2] 

• 7/kj is the electric conductance (gap junctions) between two different neu- 
rons j and k. Although, the gap junction strength is biophysically non 
static [651 I38| . we model it here as a simple ohmic conductance, tend- 
ing to equalize membrane potentials of the neurons they connect. As a 
consequence, g^y = 'gjk > 0. 

• For chemical as well as electric contacts terms of eq. (fll) the sums hold 
on all neurons, but we afford some conductances gkjXi, ^) or gkj to be 0, 
corresponding to no connection (chemical or electrical), between neuron 
j and neuron k. In this way, we can define chemical and electric network 
topologies. The mathematical results obtained in section [3] hold for any 
such network. 

2.2.2 Update of chemical synapses conductances 

Upon firing of the pre-synaptic neuron j at (discrete) time t^j"^ (w), the membrane 
conductance gkj (t) of the post-synaptic neuron k is modified as [71] : 

9k,{t) = 9kAt'f \^)) + Gk,ak,it - )(w)), t > t^pioj), (3) 

where Gkj > characterizes the maximal amplitude of the conductance during 
a post-synaptic potential. 

The function akj (called "alpha" profile [52]) mimics the time course of the 
chemical synaptic conductance upon the occurrence of the spike. We assume 
that the alpha profiles have the form: 

akj{t)=h{t)e'^oH{t), (4) 

with a(0) — 0, where h{t) is a polynomial function continuous at 0, typically 
h{t) = where Tkj is the characteristic decay time. H{t) is the Heaviside 
function. 
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As a consequence, upon the arrival of spikes at times tp (w) in the time 
interval [s,t\, eq. ([s]) becomes: 



9kj it) = gkj (s) + Gko ^ auj (t - t'f '' {uj)j . 



If we assume that this relation extends to s ^> — oo and if we set lims_i._oo gkj{s) — 
(see |T3] for a justification), we finally obtain: 



The notation (i, lo) makes explicit the dependence of the conductance upon the 
(past) firing activity of the pre-synaptic neuron j. This dependence decays 
exponentially fast thanks to the exponential decay of akj ■ 

2.2.3 Firing dynamics 

If, at time some neuron k reaches its threshold 9, Vk{t) = 9, then this neu- 
ron elicits a spike. In neurophysiology, the spike emission is a complex non- 
linear mechanism involving additional physical variables (probability of open- 
ing/closing of ionic channels), as described by the Hodgkin-Huxley equations 
[36j . The time course of the spike has a typical shape (depolarization / repolar- 
ization / refractory period) with a finite duration. One important simplification 
in Integrate-and-Fire models is to describe neuron dynamics in terms of the 
membrane potential only. As a consequence, the spike shape has to be simpli- 
fied. In the simplest Integrate-and-Fire models [311 HI] a spike is registered at 
time t whenever Vk(t) = 9, whereas the membrane potential is instantaneously 
reset to 0. Moreover, the refractory period and transmission propagation de- 
lays are set to 0. Beyond the bio-physical fact that reset, refractory period 
and transmission propagation are not instantaneous, this modelling leads to se- 
vere mathematical problems and logical inconsistencies, such as the possibility 
of having uncountably many spikes within a finite time interval, or situations 
where the state of a neuron cannot be defined (see [TH [TT] for a discussion) . 

To avoid those problems we consider that spikes emitted by a given neuron 
are separated by a minimal time scale T^ep > 0. Additionally, to conciliate 
the continuous time dynamics of membrane potentials and the discrete time 
dynamics of spikes we define the spike and reset as follows. 

1. The neuron membrane potential Vk is reset to at the next integer time 
(in S units) after t, namely [t] -\- 1. Between t and [t] + 1 the membrane 
potential keeps on evolving according to ([T]). The main reason for this 
modelling choice is that it makes the mathematical analysis simpler. 

2. A spike is registered at time [i] -1-1. This allows us to represent spike trains 
as events on a discrete time grid. It has the drawback of artificially syn- 
chronizing spikes coming from different neurons, in the deterministic case 
[141 144| . However, the presence of noise in membrane potential dynamics 
eliminates this synchronization effect. 

3. We consider that r^ep > is a multiple of 5 (thus an integer). 




r-t{''\Lj)<t 
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4. Between [t] + 1 and [t] + Tgep the membrane potential Vk is maintained to 
(refractory period). From time [t] +Tsep on, Vk evolves according to Q 
until the next spike. 

5. When the spike occurs (at time [t] + 1), the raster ui as well conductances 
gkj{t,uj) are updated. 

Figure [l] illustrates these modelling choices. 
V 




t [t]+l [t]+X, 



Sep 



Spikes 



Figure 1: Top: Continuous time course of the membrane potential in our model. 
The green dashed curve illustrates the shape of a real spike, but what we model 
is the blue curve. Bottom: A spike is registered at integer time [t] + 1. 



2.2.4 Matrix- Vector representation of subthreshold dynamics 

The sub-threshold dynamics can be rewritten in the form of a stochastic linear 
non-autonomous and non-homogeneous differential equation: 

C^+[Git,uj)-G]V^I{t,L,), (5) 

where C is a diagonal matrix which contains the capacity of each neuron. For 
simplicity we assume from now on that all neurons have the same capacity c so 
that C = cIn where In is the identity matrix of dimension N; V is the vector 
of membrane potentials; G{t,uj) is a diagonal matrix: 

G is a matrix with entries gkj for k ^ j and — gkj for k — j. It is therefore 
symmetric and the non diagonal part of this matrix specifies the connection 
topology of electric synapses in the network. Finally, I{t, u) is the vector of 
currents that can be separated in 3 components. 

/(t,w) = /(^")(i,w)-f /(""*Ht)+/(^Ht), (7) 



9L.k 



N 

E 

i=i 



5ki =^ 9k{t,uj)Ski- 



(6) 
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with: 

4=^)(f,c.)=^iyfe,afc,(t,u;), (8) 

3 

where: 

Wk,=Gk,E,, (9) 
is the synaptic weight from neuron k to neuron j] 

4'^^*)(t)=.9,,,i?,+.f*\t); (10) 

and 

The first term corresponds to the current received by neuron k from chemical 
synapses (cs) , the second term represents the external current (ext) and the last 
term (B) is the the noise part of the current. 

Define: 

^{t,u) = C-^{G-G{t,uj)), (11) 

which is a symmetric matrix. All entries of <i>(t, uj) are bounded and continuous 
in time (alpha profiles are continuous in time and the components of $(t, lo) are 
composition of continuous functions). 

Defining: 

f{t,Lo) = C-^&'\t,oj) + C-^&=''\t), (12) 

and using the decomposition of currents Q the system ^ can be expressed 
as a system of coupled Stochastic Differential Equations (SDE) of Ornstein- 
Uhlcmbeck type (0-U) in 

^^C-\G- G{t, uj)) V + C-^&') {t, Lj) + c-^ii-^^i) +^-1 (i)^ 
at v ^ v ^ / 

dV = mt,uj)V + f{t,uj))dt+ ^lNdW{t), 

(13) 

V{to) - V, 

where v is the initial condition at time tQ. Here f[t^uj) is thus a non random, 
measurable and locally bounded function of t; — is a constant; and W{t) is a 
standard A^-dimensional Brownian motion independent of i; G . 

2.2.5 Remarks 



• Although ( 13 1 is a linear system, it has a rather complex structure, due to 
the w-dependence of $(t, w), /(t, uj). Indeed, these functions integrate the 
past spike-activity of the network from the initial time to to time t. Thus, 
the membrane potentials at time t are determined by the past spikes- 
activity which, in turn, is determined by the trajectory V{s),s £ [tt),t\ 
of the membrane potentials: the evolution of the network depends on its 



whole history via w. If w is given, the integration of ( 13 ) generates a flow 



which is explicitly computed in section 3.1.4 
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• The discrete structure of the spike trains set (discrete events and discrete 
times) induces a partition on the set V of trajectories of V . A trajectory 
V & V belongs to the partition element V^i-, associated with the spike train 
w, if and only if: 

{a;fe(n) = maxig]„_i.„] Vu(t) < d; 
(14) 
ujk{n) ==1 <^ niaxtg]„_i^„] Vk{t) > 9; 

This constitutes compatibility conditions between spike trains and mem- 
brane potential trajectories. 

In the absence of noise (cr^ = 0) some partition elements V^j (depending on 
model-parameters) are not visited by any trajectory. It has been shown in 
some variants of ( jls] ) considered in [TUIIII], that the set of non empty P^iS 
is finite, leading to specific, although quite rich, periodic orbit structure 
of the attractors. In the presence of noise, all VuiS are visited by any 
trajectory with a positive probability. 

• The Wiener process on noise trajectories induces a probability measure /i 
on rasters, described in section [4j characterizing spike train statistics. 

2.3 Examples 

The model ([T]) is quite general as it considers chemical and electric synapses, 
with a spike-history dependence. We don't know any analysis of this model in 
its most general form in the literature. However, upon simplifications it reduces 
to several models studied in the literature. 



2.3.1 No electric synapse 

As mentioned in the introduction the model without electric synapse has been 
studied in [T3]. The sub-threshold dynamics reads: 

= -9L,kiVk - El) -Y,9kjit,u){Vk - E^) + h{t). 

3 

In this case 

$(t,a;) = C-^G(t,Lo), 



while f{t,uj) takes the same form (12 1 as in the general model. 



2.3.2 Simple chemical conductances 



By "simple" we mean a conductance matrix that takes the form G{t,u)) = 
K{t,uj)lN, where K{t,Lu) is a real function. Therefore, at time t and given an 
history w, all neurons have the same conductance K(t, uj), and the sub-threshold 
equation ([5| reads: 

-K{t,uj)XNV + GV + l'^'"'Ht,uj)+&'''\t) + aBdW{t). (15) 



C 



dt 



On mathematical grounds being "s imple" is a general condition ensuring 
that G{t) and G commute (see section 3.1.2). On practical grounds, classical 
examples are K{t, w) = or K{t, w) = cste. 
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Prominent illustrations are the Medvedev model [S5] and the Ostojic-Brunel- 
Hakim model [55] , 



3 Solutions of the stochastic differential equa- 
tion 

We now derive several mathematical results allowing the integration of the SDE 
( 13 1 and the consideration of firing dynamics. 

3.1 Flow in the sub-threshold regime 



We consider first the integration of (13) on a time interval [to,t\ in the sub- 
threshold regime, Vk{s) < 9, k — 1,...,N, s e [io,i]- We assume that cj is 
fixed. Note that, necessary, ui obeys LOk{n) ~ 0, k = 1, . . . , N, to < n < t {ci the 
compatibility conditions (14|). Nevertheless, we don't make any assumption on 
w before time Iq, that is, we can have any spike history prior to tg. 

3.1.1 General form of the flow 

We start by solving the associated homogeneous Cauchy problem 

V{to) = V, 

The following theorem is standard and can be found e.g. in [51. 



(16) 



Theorem 1 //$(i,a;) is a square matrix whose elements are hounded, the se- 
quence of matrices Mk{to,t,U!) defined recursively by: 

Mo{to,t,uj) =Ijv 
Mk{to,t,Lj)=lN+ f <^{s,uj)Mk-i{s,t)ds, t<ti, 

Jto 

converges uniformly in [tQ,ti]. We note: 

r(io,i,c^) lim Mkito,t,ij) (17) 
the limit function called "flow". 

3.1.2 Exponential flow 

If <i>(t,w) and ^{s,Ld) commute, Vs,t, then the flow takes the form of an expo- 
nential: 

From ( pj] ) the commutation condition reads [ G — G(i, ui) ] [ G — G{s, cj) ] = 
[ G — G{s, uj) ] [ G — G{t, uj) ] . Since we have assumed that all neurons 
have the same capacity, the commutation condition reduces to [ G — G{t, w) ] [ G — G(s, w) ] 
[G — G(s,aj)] [G — G(t,w)], i.e. G and G{t,uj) commute for all ts. Since 
G(i, w) is diagonal, the commutation condition reads gi{t, oj)Gij = Gijgj(t, lo), ^i, j. 
Therefore, here are the only cases where the commutation property holds: 
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(i) G = 0; 

(ii) G{t,uj) ~ K{t,uj)XN where K{t,uj)) is a real function. 
These cases correspond respectively to the following examples. 

No electric synapse 

In this case — — ^ G{t,u!) is a diagonal matrix. Thus, the flow takes the 

exponential form: 

I (to,t,uj) — e-'*o ^ ' — e " ■'*o ' , (18) 
which is also a diagonal matrix. 

Simple chemical conductances 

In this case w) = \G — and the flow takes an exponential form: 

r(io, i, io)^e--^ ^'0 ^-e^^(*-*«) . (19) 

It is not diagonal in the canonical basis, but it can be diagonalized by an or- 
thogonal variable change. 

3.1.3 General form for the flow 

In the general case, namely in any model taking into account simultaneously 
chemical and electric synapses with a generic form, G and G{t, uj) do not com- 



mute, and the flow (17) does not read as an exponential but as a general Dyson 
series: 

r{to,t,uj) =In + / ■■■/ A(si,a;)...A(s„,w)ds„...rfsi. 

Setting B = G^^G and A{t,uj) = —G^^G{t,oj), this equation reads: 

T{to,t,u;)^lN + J2 ••• / Y[{B + A{sk,u;))d,,---ds„^, 

where Y[ denotes the matrix product, hence ordered. Finally, the flow can be 
written: 

r(to,t,w) = In 



n=l Xi = ( B, A(si, ) "'o •'to i^^i 

X2 = {B, A(s2.i^) ) 

X„ = (i3, A(a„.a.)) 



This form is quite cumbersome when A(t,uj),B do not commute. It has never- 
theless the following property. 
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3.1.4 Exponentially bounded flow 

Definition: An exponentially bounded flow is a two parameter (to,t) family 
{r{to,t,uj)}t<to of flows such that, Vw G fi: 

1. r(io, to,uj) = Ijv and T{to, t, uj)T(t, s, u) = r(io, s, i^) whenever !i ^ 5[ s; 

2. For each v e and w e il, (to,^) ^ r(^Oi^j'^)^^ is continuous for to < ^; 

3. There is ill > and to > such that : 

I |r(s, t, w) 1 1 < Me-™(*-^) ,s<t. (21) 



Recall that a strong solution of the SDE ( 13 1 is a stochastic process for which 
the paths are right-continuous with left limits everywhere with probability one, 
adapted to the filtration generated by W{t). From [92] we have the following 
theorem: 



Theorem 2 If [11) converges to an exponentially hounded flow r(io, t, lo), ttiere 



is a unique strong solution for t> to given by: 

V{to,t,u:)^T{to,t,L,)v+ f Tis,t,u:)f{s,L,)ds+^ f T{s,t,Lj)dW{s). (22) 

Thus, given an initial condition u at a time t^ and a noise trajectory, this 
equation gives the membrane potential vector at time t by integration of the 
flow, provided maxfc^i . jv max„g[fQ Vk{u) < 9 (sub-threshold dynamics). This 



is a classical form although F has a complex structure ( 20 ) and a non trivial 
dependence in the raster history. 



Let us now show that (20) converges to an exponentially bounded flow. If 



G = then F(s,i,a;) given by (18), is diagonal and exponentially bounded. In 
this case F(s,i,a;) = cZiag(e~^ J^^'"^"''^^'*") where gu is given by Conse- 
quently, setting gi^ = inf^ „ ^ ^^(m, w), the smallest conductance value attained 
when no neuron fires ever, in the absence of gap junctions, we have: 

||F(s,t,c.)|| <e-^(*-^). (23) 
This is therefore an exponentially bounded flow. 

When G 7^ we use the following perturbational result (for details see 
corollary 2.2.3 in [33]). Set $(t,a;) = -G-iG(i,a;). Considering B = G-i(G), 
as a (not necessarily small) perturbation of ^>(t,w), we have: 

The orem 3 Let the flow T{s,t^Lj) he the exponentially hounded flow oj equation 



(23), obtained when ^(t,iu) ^ —C ^G{t,Lo). Then the flow T{s,t,LL>) of equation 



(16), when G 7^ and $(t,a;) = $(t,a;) + B obeys the inequality 
|iF(s,t,w)|| < e-^(*-^)eJ'^'ll*('--'^)-*('--")ll'*^ 

Moreover, 

||F(s,i,w) - f(s,t,w)|| < e-"-^^'-'\e^s\mr^^)~'^{r,^)\\dr _ (24) 
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Using (23) the flow associated to equation (131 satisfies the following in 



equality with ||i3|| = b: 

\\T{s,t,oj)\\ < e-^(*-«)e''(t-«) ^ ^ib~^)(t-s) (25) 



Therefore we can ensure exponentially bounded flow for equation ( 13 ) when 
b < Since G is symmetric and C — cIn the norm b — ||C~^G|j is equal to 
^ , where Ci is the largest eigenvalue of G and has the physical dimension of a 
conductance. So, theorem [s] provides a (sufficient) condition for the existence 
of a strong solution, given by: 

ai<gL. (26) 
The largest eigenvalue of G has to be smaller than the leak conductance. 

3.1.5 Remarks 

• As stated e.g. in [57] the typical electrical conductance values are of order 
1 nano-Siemens, while according e.g. to [91] the leak conductance of retinal 



ganglion cells is of order 50 micro-Siemens. Therefore, the condition (26) 
is compatible with the biophysical values of conductances in the retina. 



Looking at the series ( 20 ) one may think that the exponentially bounded 
flow is ensured whenever ^(t) have a negative spectrum. This prop- 
erty is in general not determined by the eigenvalues of <l>(i) in the non- 
autonomous case. There are examples in which the matrix $(t) have 
negative real eigenvalues Vi, but the solutions of the corresponding differ- 
ential equation grow in time. For a review and intuitive explanation see 
|41) and example 3.5 in [16j . For a more complete mathematical analysis 



see [921 . Therefore, in order to ensure a unique strong solution to ( 13 ) one 
needs to find conditions ensuring exponentially boundedness. 



The Dyson series ( 20 ) is obviously intractable on numerical and on ana- 
lytical grounds. A natural idea is to truncate this series. If gap junctions 
conductances are small compared to chemical synapses, one can reorder 
the series writing first terms containing no term G, then one term G, 
two terms, and so on. The order of truncation is fixed by the norm of 
G, which is equal to its spectral radius since this is a symmetric matrix. 
If this radius is small enough truncation at small order in G provide a 



good approximation of (20). These aspects will be discussed in a separate 
paper. 

3.2 Flow and firing regime 



The flow ( 20 ) characterizes the evolution of the membrane potential vector in 



the sub-threshold regime. Let us now extend its definition to the firing regime. 
3.2.1 Extended flow 

For a raster w, recall that t^^\u!) is the r-th time of firing of neuron j in the 
raster uj. Let T^p'>{(jj) be the p-th time in which a neuron is reset, i.e. T^^-'(w) — 
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inmj{t'^^\uj)}, and T^^H^) miiij^r {t^[\io) > T^p~^\lj)}. Let IC^p^lo) be the 
set of neurons firing at time T'-'"-' (w) i.e. 

Consider a membrane potential trajectory compatible with Thus, for an 
initial time to < T^"'"-* {u}) , Vk < 6, k = 1, . . . , N . The membrane potential vector 



V follows the evolution (22 1 until time T*-^^(w), where the membrane potential 
of neurons in /C'^^(w) is set to 0. It keeps this value during the refractory period, 
until time T''^''(w) + T ^ep. During this period the neurons ^ /C'^^^(a;) keep on 



evolving according to (22), but flow r(to,s,a;) undergoes a variation at time 
T^^^u}), where conductances are updated. The smoothness of this variation 
depends on the assumed regularity of the a function After the refractory 



period, the membrane potential of all neurons follows the evolution (22), until 
another group /C'^^^(a;) fires at time T^^-'(cj), and so on. 

To take into account the refractory periods where neurons are in the rest 
state, we introduce a diagonal matrix x = xi-^i'-^) with entries: 



Xkk{s,uj) 



0, ifs e UjrM(^),TMM +T,ep],andfc e /C^H; 

1, otherwise. 



Then we replace F in ( 22 ) by Tx- This gives 



V(to,t,uj) =r{to,t,uj)xit,L^)v + / r(s,t,cj)x(s,w)/(s,w)ds 

J to 



Oj3_ 

c 



t 

r{s,t,oj)x{s,Lo)dW{s). 

to 



To simplify the computations made in the section [4l we add the following 
modelling choice. Denote Tk{t,uj) — Tgep the last timdH before t where neu- 
ron k has been reset in the past. When to < Tk{t,u) we replace the integral 
ft„^i.s,t,uj)x{s,uj)dW{s) by jl^f^^^^^T{s,t,uj)x{s,uj) dW{s). In this way, the 
stochastic dependence upon the past is reset when the neuron's membrane po- 
tential is reset. On one hand, this modeling choice does not deeply change the 
phenomenology. On the other hand it allows us to formulate the determina- 



tion of spike train statistics as a first passage problem (section 4.1). On the 

opposite, we keep the (deterministic) integral / T[s,t^u))x{s,^j-i) f{s,L>j)ds un- 

■^*o . . . 

changed. This term contains indeed a deep effect intrinsic to gap junctions, 
explained in details in section [3. 2. 2[ 

To summarize, we have thus: 

y(to,t,w) = r(to,t,c^)x(i,^)i'+ / r(s,t,L^)x(s,t^)/(s,w)ds (27) 

OB "* 



V{a,t,uj)dW{s). 

^By construction of the model this is an integer < [t] — 1 only determined by t,Lj. 
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3.2.2 Role of gap junctions on history dependence 



In the case G = 0, F is diagonal (eq. (18l) as well as T x- The reset of the 
membrane potential has the effect of removing the dependence of Vk(t,Lo) on its 
past since Vfc([t] + r^epj is replaced by 0. As a consequence, the k component 



of eq. (27) holds, from the time Tk{t,(jj), introduced in the previous section, up 



to time t. Therefore, eq. (27) factorizes as a set of N equations pl^ : 



Vk{t,uj) = Tkk{Tk{t,uj),t,Lo)vk+ I Tkk{s,t,uj) fkis,uj)ds 



Tk(t,u:) 

t 

rkk{s,t,Lj)dWk{s), (28) 



with Tkk{'i'k{t,ijj),t^Lu) — e <= ItQ Ski^.^^jds agreement with (18 1. 

In this case, the reset has the effect of erasing the dependence of Vk on its 



past anterior to its last firing time. Note that, in equation (28), the flow is 
integrated from the time Tk{t,u}), but the total conductance defining the flow 
(the term J^^ <7fc(s, uj)ds) corresponds to an integral starting from the initial time 
to- This is because, contrarily to membrane potentials, conductances are not 
reset when a neuron fire. 

The situation is quite more subtle when electric synapses are present, G ^ 0. 
Neuron k membrane potential is still reset at firing. From this time on, its 
evolution depends on the whole vector V{t), in particular Vj{t). But Vj{t) 
depends on Vk{s)s<t via the gap junction connection. Due to this interaction 
type, the evolution of Vk depends on its past before firing, via the membrane 
potential of the other neurons. 

3.2.3 Initial conditions 



Equation (27) still depends on the initial condition V{to) = v. However, we 
are free to choose any < T^^^'iu)). Especially, we can take tQ — — oo. This 
corresponds to a situation where the neural network has started to exist in a 
distant past (longer than all characteristic relaxation times in the system) and 
we observe it after transients. This corresponds to an "asymptotic" regime 
which not necessarily stationary, if the external current depends on time. 



From the exponentially bounded flow property (21 ) F(to, t, uj)v — > as 



-oo. Therefore upon taking t^ — >■ — oo we may write (27) as 



V{t,uj)^ f r{s,t,uj)x{s,uj)f{s,uj)ds+— f r{s,t,uj)dW{s). 

J —oo ^ Jrf^(t,ij) 

The integrals are well defined thanks to the exponentially bounded flow 
property. From now on we work with the limit to — > — oo. To alleviate notation 
we remove the variable to in the expression of the membrane potential. 

3.2.4 Membrane potential decomposition 

The stochastic process V{t, w) is therefore the sum: 

V{t,u}) = V'-'^^t,u;) + V^''°""'\t,uj), (29) 
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of a deterministic part: 

J — OO 

with: ^ 

= - / r(s,t,c.)x(s,w)/(^^)(s,c.)ds, (31) 

the chemical synapses contribution to the membrane potential; 

V^--'\t,uj) = - f Tis,t,uj)x{s,u:)&-'Hs,Lo)ds, (32) 

C J-oo 

the external current + leak term contribution, and a stochastic part: 

y(""^^)(i,cj) = — /' r(s,t,cj) dW(s). (33) 

Thanks to the limit to ^ which has removed the dependence in the 
initial condition V^'^'^ (t, w) is uniquely determined by the spike history w (and the 
time dependence of the external current if any). Likewise, V^"'°'^'"'\t,uj) 

is the integral of the Wiener process with a weight T{s,t,u!) depending on the 
spike history. 

4 Spike train statistics and Gibbs distribution 

This section is devoted to the characterization of spike train statistics in the 
model. The main result establishes that spike train are distributed according 
to a Gibbs distribution. Note that we do not make any assumption on the 
stationarity of dynamics: the present formalism affords to consider as well a 
time dependent external current (stimulus). There is no explicit form of the 
potential determining the Gibbs distribution in the general case, but upon an 
assumption discussed in section |4.3[ the potential can be approached by an 
analytic form. This form, is further discussed in the section [5] with its deep 
connections, on one hand, with the Generalized Linear Model (GLM) and on 
the other with maximum entropy models. 

For the non specialized readers, we give here the main ideas behind math- 
ematics. Spike statistics is characterized by a family of transition probabilities 
giving the probability to have a spike pattern w(n) given a past spike history 
w"^^. From this set of transitions probabilities, one defines the Gibbs distri- 
bution depending on the bio-physical parameters defining the model. This is 
done in a similar way as homogeneous and positive Markov chains transition 
probabilities define the invariant distribution of the chain, although the case 
considered here is more general (we do not assume stationarity and we do not 
assume a finite memory). 

4.1 Transition probabilities 

We want to determine the probability P [w(n) | w"^^ ] to have a spiking pat- 
tern u}{n) at time n given the spike history w"^^. This can be stated as a first 
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passage problem [51 121 [SS] . 



Fix oj, n and t < n. Set: 



{t,uj)^e-V^''\t,uj), (34) 



the distance of the deterministic part of the membrane potential to the thresh- 
old. Neuron k will emit a spike at time n if there exist a time t € [n — 1, n] such 
that Vl"°'''\t,u;)^9kit,uj). 
Denote: 

al{t,u) = ^Y. rlis, t,u)ds. 

Following |85] Dubins-Schwarz' theorem can be used to change the time scale 

to write V^^"°'*'^^ (t, w) as a Brownian motion V'jf"°*'*^''(i,a;) =^ W„2i^f^^^ and the 
spiking condition at time n reads: 

This equation characterizes the first hitting time hk of the Brownian motion 
1^3.2(4) to the "boundary" 9k{t,uj). 

Denote P [/ii, . . . , /ijv] the joint law of the first hitting times of neurons 
1,...,A^. For a spiking pattern a;(n) divide the set of indices {l,...,iV} 
into two subsets S^{n,u)) = { fc g { 1, . . . , } , a;fc(n) = 1 } and S~{n,uj) — 
{fee {l,...,7V},Wfc(n) =0}. Wehave: 



V[uj{n)\ujl-^]^J> fl {/ifcG [n-l,n[} n f| {hi > n) 

keS+{n.uj) leS-{n,ui) 

"(35) 

This equation can be written in terms of an integral of the join density of 
hitting times. 

The first passage problem can be solved in simple one dimensional situations 
following a method developed by Lachal [IS] . The Laplace transform of the first 
hitting time density can be obtained as a solution of a PDE. However, we haven't 
been able to find a general form for the joint density of the TV-dimensional 
problem in our model. 

An alternative approach is to use the Fokker-Planck approach for IF mod- 
els developed by several authors [3 |48l HI El [63] . The resulting equations are 
seemingly not solvable unless using a mean-field approximation in the thermo- 
dynamic limit as done e.g. by Ostojic-Brunel-Hakim 63J in the specific case 
where chemical and electric synapses are constant and equals. The Fokker- 
Planck method is a dual approach to the first passage problem but it meets the 
same difficulties. 

In the next section and in the appendix, we establish general conditions 
ensuring nevertheless that this system of transition probabilities define a Gibbs 



distribution. Then, in section 4.3 we propose an approximation allowing an 



explicit computation of transition probabilities. 
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4.2 The Gibbs distribution 



The transition probabilities (35) define a stochastic process on the set of rasters, 
where the probabihty of having a spiking pattern uj{n) at time n depends on 
an infinite past w"^^. This is the extension of the concept of Markov chain 
to processes having an infinite memory. Such processes are called "chains with 
complete connections" [S51[S11[5S]. Under suitable conditions, exposed in the 
appendix, a sequence of transition probabilities defines a unique probability 
distribution ^ on the set of rasters, called a "Gibbs distributioi^' . We give 
here its main properties. 

4.2.1 Gibbs potential 

The function: 

^{n,uj) = \ogV[uj{n)\oj^-J], (36) 

is called a Gibbs potential. Then, the probability of any block uj^^ n > m, given 
the past w™^^ is given by: 

P[< = (37) 

This form reminds Gibbs distributions on spin lattices, with lattice translation- 
invariant probability distributions given specific boundary conditions. Given a 
spin-potential of spatial-range n, the probability of a spin block depends upon 
the state of the spin block, as well as spins states in a neighborhood of that 
block. The conditional probability of this block given a fixed neighborhood is 
the exponential of the energy characterizing physical interactions, within the 



block, as well as interactions with the boundaries. In (37), spins are replaced 



by spiking patterns; space is replaced by time. Spatial boundary conditions are 



here replaced by the dependence upon the past ui'^^o ■ The potential (36) has 
important properties: 



It corresponds to a process with memory, as opposite e.g. to other sta- 
tistical physics approaches of spike trains statistics based on maximum 
entropy models of Ising type [7SJ UHl HSj and its extensions [IHl US] • In 
those models, spike dynamics has no memory and the probability of events 
occurring at different times is simply the product of probability of those 
events. 

It has infinite range, corresponding to an infinite memory. This is due 
on one hand to the synaptic response s ([4|) and on the other hand to the 



presence of gap junctions (see section 3.2.2 1. This can be readily seen on 
the form of the flow ( [20| . Fortunately, the exponentially bounded flow 
property allows to approximate the potential by a finite range potential 



corresponding to a Markov process (section 5.3.1) 



It is normalized, since is the log of a transition probability. Hence, there 
is here no need to compute the "partition function". 



^Here, Gibbs distributions are considered in a general setting. In the language of statistical 
physics we consider infinite-range potentials and non translation-invariant potentials I75II31I . 
However, since we are considering time evolving processes we shall use the (equivalent) lan- 
guage of stochastic processes and chains with complete connections |25. .52, . 



18 



• It depends on the parameters defining the network (synaptic weights, ex- 
ternal currents, biophysical parameters) in a nonlinear way, as opposite 
to, e.g. maximal entropy models of Ising type [76l |83] and extensions 
[28l [29] . where the potential depends linearly on parameters, but where 
the number of parameters grows fast with the number of neurons, and 
whose interpretation is controversial. 



4.2.2 Gibbs distribution 



In this context, a Gibbs distribution is a probability measure fj, on the set of 
rasters such that, for any block o;^: 



=Er=™ </>( '.'^ ) 



li{dio). 



(38) 



where the notation has been introduced in section 2.1 

The Gibbs distribution has the following properties (see [88l [58] for more 
details). 



When the stochastic process defined by the transition probabilities (35) 



is time-translation invariant (this is the case if the external current is 
time-independent), the Gibbs distribution satisfies a maximum entropy 
principle : it maximizes the sum "entropy + average value of the poten- 
tial" . This is equivalent to free energy minimization in statistical physics. 



This aspect is developed in section 5.3.2 



• In the case of time-translation invariant transition probabilities depending 
upon a finite past (homogeneous Markov chain) the Gibbs distribution is 
the invariant measure of the chain (it is unique provided that transition 
probabilities are strictly positive). This aspect is developed in section 

mi 

• Gibbs distributions are also defined in the case where the potential is 
not time-translation invariant and has infinite range, upon mathematical 
assumptions on transition probabilities stated in the appendix. 



4.3 Approximating the Gibbs potential 

The previous subsections provide mathematical results, but they do not allow 
explicit computations. The basic step for this would be to have an explicit form 
for the Gibbs potential (resp. the transition probabilities). As mentioned in 
section [4.1| this is not tractable in general. We propose here an approximation. 



4.3.1 The stochastic term y("°"^)(t,aj) and compatibility conditions 



y (noise) ^-j jg ^Yie integral of a Wiener process (eq. (33 1). As a consequence 
this is a Gaussian process, with independent increments, mean 0, and covariance 
matrix: 



Q{t, uj) =^ Cov 
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where denotes the transpose. From standard Wiener integration we have: 



r{s,t,u;)x{s,io)dW{s)i / r{s',t,uj)x{s,uj)dW{s')) 

Tk(t,U}) JTk(t.,Uj) 



„2 

^ JTk(t,U}) JTk{t,Ul) 

SO that 



r(s, t, io) xi-s, uj)x^is', c^)r^(s', t,uj)E[ dW{s)dWi.s') 



where we used x{St^)x^ {s,uj) — x(s,i^)- 



(39) 



However, as mentioned in section |2.2.5| the model definition imposes com- 
patibihty conditions between the trajectories of the membrane potential V{t, u) 
and the raster a; (see eq. (14)). Fixing a raster fixes the flow F as well as 
the deterministic part of the membrane potential V'^'^'^ {t^uj) (eq. (30)). Then, 
y(nozse)^^^ has to obcy the following compatibility conditions: 



Wfc(n)=0 <^ max„gi„_i „i Vl"'"'""' {u,u) - ek{u,u) 
yk^l,...N,yn<t, { ■ . ^ . 

ujk{n) = l ^ max„g]„_i^„] V^"'"''"'' {u,uj) - dkiu^uj) 

(40) 

As a consequence, when computing the law of V^f"°*'"^-*(.,a;) we have to 



take these constraints into accoun 



of K 

Pr 



(noise) 



We write Pr 



(noise) 



(noise) / 



the law 



(.,w) given these constraints. Although, V'l'^) is Gaussian, 

^(nojse)^^^^ is not a Gaussian distribution. For example, if U!k{n) — 

then, necessarily, V^"'°^^'^\t^uj) < 9 ~ V^'^\t,LL!), [t] + 1 = n, while a Gaussian 
random variable takes unbounded values. 



4.3.2 A Gaussian approximation. 

A plausible approximation consists then of approximating Pc [ T^*^"°*^'^-'(., oj) ] 
by the Gaussian law of ^("o*^^)^ j "ignoring" compatibility conditions. This 
approximation can be justified under the following conditions. 



• Weak noise. Consider the compatibility conditions for Wfc(n) = 0: Vu £ 
]n — Vj|"°*''^^(u, w) < 9k{u,uj). What is the probability that the 

Gaussian noise V^^°^^'^\u,uj) violates this condition for some u in the 
interval ]n — l,n\ ? 

Denoting cr^(u,a;) Qkk{u,uj), the variance of V^f"""'^'' (u, w) at time u 
(given explicitly by eq. (39)), this probability is given by: 



P 



(noise) 



{u,uj) > 9kiu,uj) 



+ 00 



fc(".-) 



-^dh. 



^We warmly acknowledge E. Shea-Brown for this deep remark. 
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From the explicit form of (t^(u,w) (eq. (39l) and from the exponentially 
bounded flow property of F, af.{u,u}) is upper bounded, uniformly in u and 

2 

by -^B where B is some constant depending on the model parameters 
and where as is the Wiener noise intensity. Thus, cr^(u, w) — > as ctb — > 0. 

If (Jb = then the compatibility conditions for Wfe(n) = is Vm G]ri — 1, n], 
9k{u,uj) > (i.e. V^'^\u,uj) < 9). Assume that 6k{u,uj) > e for some 
e > 0. 

If UB ^ 0, > so that: 



P 



< I e dh^ -\ I - erf , ^ 



where erf{x) admits the following series expansion as x — > +oo: 
1/1 1 3 15 



ev f (^x^ - 1 — 6 I 

•v/tt V X 2x^ Sx'^ 



Therefore, as ctb — 0, we have 
P 



l_^2j_/l_^ 3 15 



(41) 



with a; = — . 

This shows that the probability that the noise violates the compatibil- 
ity condition decreases exponentially fast as ctb — 0. As a consequence. 



when the noise is small, approximating 



by the Gaus- 



sian law of provides a reliable approximation. This amounts to 

considering that the spikes arising between [n — l,n[ are determined by 
the deterministic part of the membrane potential, not by the noise. If a 
neuron is about to fire at time n because its deterministic part crosses the 
threshold at a time t G [n — 1, n[, the (weak) noise can affect the time t 
where the crossing occurs, but, with high probability this time stays in 
the interval [n — 1, n[. 

Time discretization. Beyond the compatibility conditions, an addi- 
tional aspect makes the analytic computation of transition probabilities 
delicate: membrane potential time evolution is continuous. Since we are 
focusing here on spike statistics, where spike occur on discrete times, this 
obstacle can be raised using the following approximation. We replace the 
spike condition of neuron k at time n: 3u e]n — 1, n], V^^"°''"^^(u, a;) > 
9k{u,ui) by: V^^""'*"^"* (n, oj) > 9k{n,ui). The argument supporting this 
choice is the following. Suppose that vl."'°''^^\n — l,tj) < 9k{n ~ if 

vl^''^''^\n,u;) > 9k{n,uj) then with probability 1, the continuous process 
crosses the threshold between Jti — 1, n], therefore the spike is registered at 
time 71. On the other hand if y^("°"^) (^ji^ < 0^.(^n, uj) there is nevertheless 
a positive probability that the continuous process crosses the threshold at 
some time t between ]n — 1, n[. In this case, there will be a spike that our 
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approximation neglects since V^"'°^^'^\n,uj) < 9k{n,uj). The probability of 
occurrence of such an event can be explicitly computed for general diffu- 
sion processes (see [64]) and depend on the noise intensity as and time 
step discretization 5. It is given by: 



Cs-exp 



-2{9k{n -l,uj)- l,c.))(^fc(n,c.) - vl^'''''\n,oj)) 

San 



where Cg- is a bounded function depending on the value of 9 at time n — 1. 
Therefore, when aB,S are small (as in our case) this probability is also 
small. 

4.3.3 Approached form of the potential. 

Thanks to these approximation we have now the following result. If we note: 

]-co,9k{n-l,uj)l ifwfe(n) = 0; 

[9k{n-l,uj),+co[, ifwfe(n) = l; 



and 

N 



Jin,u;) = llJkin,Lj), (42) 



k=l 



where Y[ denotes the Cartesian product of intervals, eq. (35) becomes, using 
the Gaussian approximation for y("°''*'^)(.,a;): 



P [a;(n) | u^^-J ] = / — ^— — ^d., (43) 

Jj(n,u) (27r) 2 |Q(n- l,a;)|2 



where dv = Yl^^i dVi. 



Taking the log of ( 43 ) we obtain the approached form of the Gibbs potential. 



4.3.4 Remarks 

The apparently simple form (|43|) hides a real complexity. 



• The integration domain J^{n, uj) corresponds to a product of intervals in- 
volving the variables 9k{n — l,w), the distance of V^'^\n — l,w) to the 



threshold. Now^from eq. ( 13 ) , ( 30 ) , involves the chemical synapse current 



/'^'^''^ (i, cj), eq. (jsj), and the external current 

/(^"=*)(*) (M' integrating, via 
the flow r, the spike history of the network. As a consequence, the transi- 
tions probabilities depend on all parameters defining the system: synaptic 
weights Wkj (eq. (|9|, gap junctions, external current, and biophysical 
parameters such as membrane capacity, or characteristic time scale of the 
post-synaptic potential akj{t) (eq. Q). 

• Without electric synapses the covariance ([39]) takes a diagonal form since T 
and X S'l'fi diagonal matrices. Thus, in this case, neurons are (conditionally 
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upon the past), independent. As a consequence the transition probability 
P [w(ri) I ] factorizes: 



N 



P [a;(n) \u;Tj]^Y[¥ [co.in) \ c.Tj ] . (44) 

k=l 

5 Consequences 

In this section we adopt the following point of view. Assuming that the model 
presented here captures enough of the biophysics of real neurons, what are the 
-relevant for neuroscience- consequences of the mathematical results developed 
in the previous sections ? We essentially focus on spike trains analysis and argue 
that: 

1. Spikes correlations are not only due to shared stimulus: there are corre- 
lations induced by dynamics, that persists without stimulus. Moreover, 
in the absence of electric synapses neurons are conditionally independent 
upon the past. 



2. The Gibbs potential ( 36 ) or even its Gaussian approximation ( 43 1 , includes 
existing models for spike trains statistics such as maximum entropy models 
[761 [841 [801 mi EH IMlEa [88] and GLM [HI [Ml [B [ZQl [51] . 

5.1 Correlations structure 

5.1.1 Transition probabilities do not factorize in general 

This can be illustrated in the Gaussian approximation. The covariance matrix 
Q{n— 1,0;) can always be diagonalized by an orthogonal variable change P{n — 
1,0;) depending on n and o;: 

Q{n- l,uj)P{n-l,uj) = Pin - l,o;)E(n - l,o;), 

where Yi{n — is diagonal with real, positive eigenvalues <y\{n — l,o;). 

Upon the variable change V = P{n—1, uj)V' , the transition probability reads 
therefore: 



P [o;(n) \u;--J]^ [ f[ ^^I^ dVl 



where J'{n,uj) is the image of J{n,Ljj) in the variable change. As a classi- 
cal result the variables change has transformed the join Gaussian density of 
y{noise)(^^ — l,o;) iuto & product of one dimensional Gaussian with mean zero 
and variance (T^(n — l,o;). However, the same variable change has transformed 



the domain J'{n,u)), reading as a product of intervals (eq. 42) into the do- 
main J''{n,u)) which is not a product any more, except in the case with no gap 
junctions. 

Therefore, in general: 



N 



P [u;{n) \c."_-J]^l[F [c.,{n) \ c."_-J ] . (45) 
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5.2 Correlations 



Note first that conditional independence upon the past does not mean indepen- 
dence. Conditional independence upon the past means that the Gibbs potential 
reads: 



N 



k=l 

This is a sum of per-neuron potentials 4'n,k, but each of this potential is a 
function of the past network activity uj. On the opposite, neurons independence 
would mean that: 

N 

0n (w) = ^ (l)n,k (uJk) , 
fc=l 

where now the per-neuron potentials (pn^k depends on the past activity of neuron 
k only. 

Now, what are the sources of general correlations (i.e. pairwise and higher 
order) in our model ? Namely, what makes the potential depend on the network 
history ? Even in the absence of electric synapses, correlations can be induced 
on one hand by the stimulus (the term V'^^*^{t,uj), eq ( |32| )) if the external cur- 
rent i^'^^^\t) in (1) has correlations between its entries ^^^^\t) and i'f'^^\t) for 
k 7^ I are correlated), and on the other hand by the chemical synapses term 



V^'"'\t,uj), (eq. (311). Even if the stimulus is zero, the chemical synapses term 
remains. We arrive then to the (somewhat obvious) conclusion that in the model 
the main source of correlations is not a (shared) stimulus. It is due to interac- 
tions between neurons. 



A deeper question is to quantify the intensity of correlations induced by (i) 
chemical synapses; (ii) electric synapses; (iii) stimulus and under which con- 
ditions (parameters value) shared stimulus correlations are dominant. This is 
the main additional issue to investigate neuronal encoding by a population of 
neurons. This requires an extended investigation beyond the scope of this paper. 



5.3 The Gibbs potential form includes existing models for 
spike trains statistics 

5.3.1 Markovian approximation 



The exponentially bounded flow property ( |21[ ) means that the norm of the flow 
r(s,i, w) decay exponentially fast as t ~ s grows. A consequence of this is the 
exponential decay oj the variation ofT{s,t,uj) with respect to uj. This property, 
mathematically defined in the appendix, means in words the following. Assume 
that we know the spike patterns of the raster w from some integer time m < t 
to time t, but we ignore the spike patterns occurring before m. What is the 
maximal error than we can make on the value of T{s,t,u)) ? As shown in the 
appendix, thanks to the exponential bounded flow property, this error decays 
exponentially fast as t — to grows. The same property holds for the membrane 
potential. It holds as well for the Gibbs potential under conditions stated in the 
appendix. 

As a consequence the norm of the difference between the exact Gibbs po- 



tential ( 36 ) and an approximate potential where the past spike history 



n-l 
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is replaced by w" ]^ for some integer D > called memory depth, decreases 
exponentially fast with D. 



Therefore, we may replace the infinite range potential (36) corresponding to 
a process with infinite memory, by a truncated potential, corresponding to a 
Markov chain with memory depth D. In this case, the Gibbs distribution is the 
invariant probability of the corresponding Markov chain. For a small number of 
neurons this probability can be characterized using transfer matrices techniques 
[55] , For larger networks Monte Carlo methods can be used [55] , 

5.3.2 The maximum entropy principle 

Assume that transition probabilities (|35|) are time-translation invariant, i.e. the 



Gibbs potential (36 1 obeys 4>{n + k,uj) = (f){n,Lo), V/c. This corresponds to the 
physical concept of stationarity. Then obeys a maximum entropy principle. 
There is a function V [((>], called topological pressure, which is analogous of the 
free energy density (especially, this a generating function for cumulants) which 
obeys 

r[^]= sup {h[iy] + ly [<!>)]) ^h[f,] + (46) 

where h is the entropy rate (or Kolmogorov-Sinai entropy), Mmv is the set of 
time-translation invariant probability measures, and fi is the Gibbs distribution. 

If (p has finite range D: (l){n,uj) = (t>{io^_]j) then, a classical result from 
Hammersley and Cliff'ord [35] establishes that <j) can be written as: 

fc=0 

where the /3kS are real parameters, (non linear) functions of the network param- 
eters of the model. Moreover, the OkS are functions of type: 

Okioj''i.D) ^ ^^lih) ■ . ■OJ^„itn), k G { 1 , . . . , iV } , £ {n-D,...,n}. (48) 

These are characteristic functions of events: Ok{tUn-D) — 1 whenever neuron ii 
fires at time ti, . . . , neuron ii fires at time ti in the raster cj. It is easy to show 
that there are such functions. Note that this result is general and does 
not depend on the Gaussian assumption of section \4-.3.S\ Only the dependence 
of the PkS in network parameters is constrained by the analytical form of the 
potential. 



Using (471, eq. (46) reads: 



k=0 / k=0 



Assume now that we are only given a raster w generated by our network, and 
that we want to approximate the exact Gibbs distribution fi from those data. A 
classical approach, dating back to Jaynes [40 consists of fixing a certain set of 
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observables and to measure their empirical average value. From Hammersley- 
ClifFord theorem any observable can be written as a linear combination of ob- 
servables Ok, so there is no loss of generality in assuming that observables are 



functions of the type (48). Call Ak the empirical average of Ok- The AkS are 



fixed from the experimental data and therefore constitute constraints to our 



problem. Therefore, the maximal entropy principle (46) reads: 



2 J V 2 

P[0]= sup h[i^] + ^l3kAk = h[fi] + ^/3kAk. 

The Gibbs distribution /i maximizes the entropy rate under the constraints 
fj, [Ok ] = Ak- These constraints fix the value of the /3feS in a unique way. 

Assume now that we are not able to measure the empirical average of all 



OfcS: as the number of events implied in the product (48 ) increases its probability 
becomes more and more difficult to estimate. So, typically, one has only a subset 
of constraints. The application of Jaynes' strategy consists of maximising the 
entropy rate under this subset of constraints. The Gibbs distribution obtained 
this way is only an approximation of fi. 

The Jaynes strategy has been applied by many authors to the analysis of 
spike train statistics, especially in the retina [76 l l84 l l80 t [8T | [6T | [28 l [2 ^ [88]. 
In particular, fixing firing rates and the probability of pairwise coincidences of 
spikes leads to a Gibbs distribution having the same form as the Ising model. 
This idea has been introduced by Schneidman et al in [75] . They reproduce accu- 
rately the probability of spatial spiking pattern. Since then, their approach has 
known a great success (see e.g. [80l [TH [82] ) , although some authors raised solid 
objections on this model [72l[8Tl[84l[6T] while several papers have pointed out the 
importance of temporal patterns of activity at the network level [5ni[Sni[ZS|- As a 
consequence, a few authors |531I39[[75] have attempted to define time-dependent 
models of Gibbs distributions where constraints include time-dependent corre- 
lations between pairs, triplets and so on [SH]. As a matter of fact, the analysis 
of the data of [76^ with such models describes more accurately the statistics of 
spatio-temporal spike patterns ^88j . 

Now, the inspection of the Gibbs potential in our model leads to several 
strong conclusions: 

1. The Gibbs potential is definitely not Ising, and is actually quite far from 
an Ising model. The main reason for that is that Ising model involves only 
instantaneous (pairwise) events, with no memory effects. On the opposite, 
strong memory effects exist in our model, requiring to consider transition 
probabilities depending on spike history. 



2. The exact expansion (47) involves 2^^ constraints, making rapidly in- 
tractable any numerical methods attempting to match all constraints. 
Obviously, one can argue that some /3kS can be close to so that the 
corresponding constraints can be ignored in the approximation of /i. Un- 
fortunately, this argument does not tell us which are the negligible terms. 
It might well be that the answer depend on the network parameters as 
well. This aspect is under current investigations. 
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3. Although the expansion (47) of a potential into a linear combination of 
observables is quite natural in the realm of statistical physics and Jaynes 
principle, it may not be appropriate for the study of neural networks 
of the type studied here. Indeed, the parameters /3feS, whose number 
increases exponentially fast with the number of neurons, are redundant: 
they are functions of the network parameters, such as synaptic weights, 
whose number increases, at most, like iV^. Using approximations like 
the Gaussian one allows to obtain an explicit form for the potential as a 
function of these parameters. One obtains in the end a nonlinear function, 
but with quite a bit less parameters to tune; additionally those parameters 
have a straightforward biophysical interpretation. This constitutes an 
alternative strategy to estimate spike statistics from empirical data. 



5.3.3 The GLM model 

The Linear-Nonlinear (LN) models and Generalized Linear Models (GLM) 
[551 [551 1551 [711 [Ml m [TO] constitute prominent alternatives to MaxEnt models. 
We focus here on the GLM. 

Let X = x{t) be a time-dependent stimulus. In response to x the network 
emits a spike train response r. This response does not only depend on x, but 
also on the network history H. The GLM (and LN) assimilate the response r 
as an inhomogeneous Poisson process: the probability that neuron k emits a 
spike between t and t-\-dt is given by Afc(t) dt, where Xk(t) is called "conditional 
intensity" . In the GLM this function is given by: 

Xk{t) = f (^bk + {Kk*x)(t) + Y,iHkj*r,)(t)^ , (49) 

where / is a non linear function (an exponential or a sigmoid); bk is some con- 
stant fixing the baseline firing rate of neuron k; Kk is a causal, time-translation 
invariant, linear convolution kernel that mimics a linear receptive field of neuron 
k; * is the convolution product; Hkj describes possible excitatory or inhibitory 
post spike effects of the j th observed neuron on the k th. As such, it depends 
on the past spikes, hence on uj. The diagonal components H^k describe the 
post spike feedback of the neuron to itself, and can account for refractoriness, 
adaptation and burstiness depending on their shape. 



The GLM postulates that, given the history "H and stimulus x, the neurons 
are independent {conditional independence upon the past and stimulus). The 
network response at time n is a spiking pattern Lo{n) while the history is the spike 
activity w"^^. As a consequence of the conditional independence assumption 
the probability of a spiking pattern given the history of the network reads: 

N 

P[c.(n) \u"_-J]^l[\k{nr''^''\l-Xk{n))'--''^"\ (50) 

k=l 

and the Gibbs potential is: 

N 

(/)„(a;) = ^[wfe(7i) logAfc(n) + (l-Wfe(7i))log(l-Afc(n))]. (51) 

k=l 
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It is normalized by definition. 



Although the form of Xk{n) is postulated, the GLM has been proved quite 
efficient to analyse retina spike trains. In our model, one can have an explicit 
form for A/j(n) in the case without gap junctions and under the Gaussian ap- 
proximation. It has the form (|12j): 



Xkiuj) = 7r{Xk{n - 1,0;)), 



(52) 



where 



tt{x) = 

is the sigmoid function, whereas 



1 



2tt 
Xkin 



l,w) 



gfc(n-l,i^) 
a- ^(ra- l,i^) ■ 



There is a clear 



similarity between equation (52) and (51) where / in (49), usually an exponen- 
tial is replaced by tt, a sigmoid. In our model the stimulus and the history are 
contained in the term Xk{n — l,a;) respectively via V'^'^^^^n — l,a;) (eq. (32 1), 



the external current contribution, and V^'^'^^n — l,io) (eq 
synapses contribution to the membrane potential. Therefore 
Kk * X and Y^. Hkj * r 



31 1), the chemical 
here the terms 



in eq. (49) are made explicit. 



Note that we have 

nowhere specified the explicit form of the external current it can have 

the form {Kk * x){t) including, e.g. the integrated effect of an input x, filtered 
by neurons layers before arriving to our model, modelling a final layer. The 



term J2j Hkj 
interactions. 



corresponds, in our case, to the effect of chemical synapses 



As stated several times in this paper, the conditional independence assump- 
tion does not hold any more in the presence of gap junctions. So, in our model 
with gap junctions, the GLM is an incorrect model. 



6 Conclusion 

In this work we have analyzed the join effects of spike history dependent chemical 
synapses conductances and gap junctions in an Integrate and Fire model. We 
have pointed out several technical and conceptual difficulties mainly lying in 
the fact that, in general, the chemical conductance matrix G{t,uj) and the gap 
junctions matrix G do not commute. As we showed, this has no impact on the 
well posedness of the model, provided that the values of chemical and electrical 
conductances are compatible with biophysics. There exists a strong solution of 
the stochastic differential equations. However, except in some specific cases, the 
fiow reads as a Dyson series, hardly tractable. We discuss some approximations 
below. 

We have also considered the statistics of spike trains in this model and 
showed that it is characterized by Gibbs distribution, time-dependent (non sta- 
tionary) whenever the external current is time-dependent. The corresponding 
potential can be computed under a Gaussian approximation: it has infinite 
range with exponentially decaying interactions. 

*Note that the exponential form is used for its convexity property, useful to learn the 
parameters of the GLM [l] . 
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The main observation resulting from our analysis is that spike statistics is in- 
decomposable. The probability of spike patterns does not factorize as a product 
of marginal, per-neuron, distributions. This effect is enhanced by the presence 
of gap junctions. As a consequence, in that model, there is absolutely now way 
to defend that neurons act as independent sources. Additionally, correlations 
mainly result from the chemical and electrical interactions between neurons 
(correlations persist even if there is no external current / stimulus). Our work 
suggests especially that electric synapses could have a strong influence in spike 
train statistics of biological neural systems, especially the retina where gap junc- 
tions connections are ubiquitous. 

There are ongoing debates on the correlations between ganglion cells in the 
retina: are ganglion cells independent encoders or do they act in a correlated 
way [Sg [Ml Ea El Eg EH |29|. What are the sources of correlations [67l |77| ? 
Are they only due to a shared stimulus or shared noise, or do they also result 
from interactions between neurons |87| ? 

Our model involves only firing cells while most cell type in the retina (amacrine, 
horizontal, bipolar) do not "spike" . Nevertheless, the conclusion of our paper 
raises questions about retina spike trains: if dynamical correlations are so im- 
portant in this simple model, how can they become weaker in a more complex 
neural system like the retina ? If such correlations exist why are they so difficult 
to exhibit in experiment, without controversy |76l 172] ? Concerning this last 
question, note that our spike trains are not binned. Binning data at a charac- 
teristic time scale of 10 — 20 ms, i.e. larger than the characteristic time scale of 
gap junctions and of order the time decay of chemical synapses, could have a 
strong impact on statistical inference. The effect of binning in our model needs 
however to be further investigated. 

We have also pointed out that the Gibbs potential obtained in our model 
is quite more complex than Ising model or similar models used in retina spike 
train analysis [76l HH [28l [29l [87] . Especially, it involves spatio-temporal spike 
patterns. Handling spatio-temporal events in Gibbs distributions requires more 
subtle algorithms than simultaneous events as those considered in the Ising 
model [7B1[TS1[53] and its extensions ^H1[2S]- These aspects have been developed 

in [HH1[SH1- 

As mentioned in the introduction, one of the main motivation of this work 
was to better understand how ganglion cells in the retina supporting both chem- 
ical and electric synapses coordinate spatio-temporal spike patterns to encode 
information conveyed to the brain. Developing a further understanding of the 
regulation of gap junctions, as well as the dynamic relationship between elec- 
trical and chemical transmission, is an important challenge for the future @]. 
Although our results are still largely lying on a theoretical ground, we hope that 
this paper will help the neuroscience community to make a step forward in this 
direction. 
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7 Appendix: Existence and uniqueness of a Gibbs 
distribution 



We give here the mathematical definition and proof of existence of a Gibbs dis- 
tribution in our model. We also comment under which conditions the uniqueness 
could be proved. For mathematical references on Gibbs distributions and chains 
with complete connections see e.g. [5T1 [5^ . 

7.1 Definitions 

For n e S, we note -4"^^ the set of sequences w"^^ and J^<„_i the related a- 
algebra, while J- is the cr-algebra related with ft. J^) is the set of probability 
measures on {fl,J-). 

System of transition probabilities 

A system of transition probabilities is a family { P„ }„gs of functions P [ | ] : 
A X A'l^ [0, 1], such that for every n £ Wj, 

• for every uj(n) G ^ , P [uj{n) | . ] is measurable with respect to J-<n-i] 

. for every u:'^-J e AT^, E.in)eA^ H^i) \ uj^J ] = 1. 

Chain with complete connections 

A chain with complete connections is a probability measure fi in V(yi,T) such 
that for all n G S and all J^<„-measurable functions /: 




Continuity with respect to a raster 

We note, for n e S, m > 0, and r integer: 

Lu Lu' if = Vj' G { n — m, . . . , n } . 

A function f{t^Lj) is continuous with respect to lj if its m-variation: 

var^[f{t,.)] =sup| \f{t,uj)-f{t,uj')\ :a;'°=' w'j. 

tends to as TO — >■ +00. 

Gibbs measure 

A probability measure fi in V{fl,T) is a Gibbs measure if: 

• /i is a chain with complete connection with respect to the system { P [ a;(n) | uj 
. Vn G yuj^l^ G Al-J, P [c.(n) | Lo^rJ ] > 0. 

• The system { P [uj{n) | w"^^ ] } is continuous with respect ui. 
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7.2 Existence of a Gibbs measure 



This is established by a standard resuh in the frame of chains with complete 
connections stating that a system of continuous transition probabilities on a 
compact space has at least one probability measure consistent with it [55]. We 
give therefore here the proof of the continuity of P | w"^^ ] with respect 

to w. To prove this it is sufficient to show that P | w"^^ ] depends con- 

tinuously on the only functions V^'^^n — l,a;) and Q{n — and that these 
functions are continuous with respect to w. 

7.2.1 P [w(n) I w"^^ ] depends continuously on V'^'^\n—\,u}) and Q(n — 

This property is straightforward in the Gaussian approximation. For more gen- 
eral diffusion processes the argument is based on the differentiability, and, con- 
sequently, on the continuity of the first hitting time density with respect to time 
dependent boundary perturbations |23j . This result holds for multidimensional 
correlated processes, under conditions ensuring uniqueness of strong solutions 
(exponentially bounded flow in our case). In our case, the boundary depends 



continuously on y('^)(t,w) (eq. 34) thus, the continuity of the first hitting time 



density and therefore of the transition probabilities (eq. 35 1 with respect to 
boundary perturbations implies the continuity of the transition probabilities 
P [w(n) I w"^^ ] with respect to if V^'^^{1i,,uS) is continuous with respect to w. 
Additionally, in our case, the covariance Q(t, a;) of the process 

y(no»se) depends 
also on io so we have to check that Qit^iS) is continuous w.r.t w. 

7.2.2 Useful notations and results 

In the next section we prove the continuity of V'^'^^it.iS) and Q(i,a;) w.r.t 

For this, we introduce the following (Hardy) notation: if a function /(t) is 

bounded from above, as i — -l-oo, by a function g(£) we write: /(f) -< g{t). 
We use the following bound obtained in [13]. If the alpha profile has the 

form Q where h(t) is a polynomial of degree d, then: 

Tn — 

varm[akj{s, .)] ^ 2Pd( — )e , (53) 
where Pd is a polynomial of degree d, and: 

N 

var^[gk{s,.)] ^ 2 V GfciPd( — )e"^ , (54) 

— Tkl 

1=1 

uniformly in s. 

7.2.3 Continuity of V^f*^*^ v^rith respect to uj 



Fix t, ■m>0 and lu,u!' with uj cu' . From (|9|), (|31|), we have: 



Vt\t,u:)-Vt\t.u.')\^ 
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1 /■* ^ r 



ds 



1 /■ 

1 ^ /■* 



We have: 



+ sup[rfcj(s,t,w)xjj(s,a;)] t;ar™[aj,,-'(s, .)], 

and 

var„i[Tk]{s,t, .)xj]{s, .)] < war„i[rfcj(s, .)] sup Xj]{s,^^)+sirpTkj{s,t,uj) varm[Xjj{s,i^)]- 



From equation (24 1 with $(t,w) = B + A{t,uj) and since — \G{t,uj) 



is diagonal we obtain the fohowing bound: 

\TuAs.t,uj)-Tu,{s,t,uj')\ < \\Tis,t,co)-Tis,t,Lo')\\ < e(''-^)(*-^) 
with: 



J* \\G{r,u^)-G{r,u;')\\dr _ ^' 



\\G{r,uj) - G{r,uj')\\ < N maxvar^[gk{r, .)] =^ Nvar^[gf^{r, .)]■ 

k 

Therefore, 
We have; 



t N ^ 

var^[gf^{r, .)]dr ^ 2 (i - s) ^ G^,Prf( — )e"^ , 



which converges to as to +cxi, so, that: 



c ^ki 



N 



Therefore, 



N 



1=1 



N 



[rfe,(s,i,.)]^2-(t-s)e(''-^)(*-«' ^G^jP,( — )e . (55) 

1=1 



As shown in |13| ajji{s,uj) is bounded Vj, j' by a constant a+ < oo while, 
following the definition of x(s,u;) given in section 3.2 sup^g^j xl*;"^) = 1- 
Moreover, sup^g^ Tfcj (s, w) < sup^g^ \\r{s,t,u)\\. Therefore: 

var„^[rkj{s,t, ■)xjj{s, ■)ajj'{s, .)] 
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< 



varjn[Tkjis,t, .)] + sup \\r{s,t,u})\\varm[Xjjis, ■)] 



a++sup \\r{s,t,uj \\varm[ajj'{s, .)], 
wen 



and applying inequalities ( 25 ) and i53l we obtain: 



var^[Tkj{s,t, .)xjjis, ■)ajj'{s, .)] ^ 



Therefore: 



N 



N 



1=1 ki 



or / TIT- \ — — 

^ 'jj' 



e(''-^)(*-)d,+ 

) 



2P,( — )e 
'''jj' 

N 



1 = 1 



a 



varm[Xj]{s, .)]ds] 



The integrals are convergent when 6 — — < 0, i.e. when the condition ( 26 ) 
is fulfilled. In this case: 

1 ^ 



c(6 



(=1 



e 'fc! + a 

'kl 



varjn [Xjjis, ■)]ds 



Taking m — +oo the two first terms converge to exponentially fast. Con- 
sidering the third term we have: 



var„i[Xjjis, ■)]ds 



The second term is zero since varjnixjji^^ .)] = on this interval, by definition. 
For the first term, varm[xjj{s, ■)] < 2 so that: 



e^''-^'>(*''')var^[Xni^,-)]ds< 
which converges to exponentially fast as m — )■ 



c 

f oo if 6 — ^ < (condition 



Therefore, V^'''^ is continuous with respect to w with an exponentially fast 
decay of its variation. 
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7.2.4 Continuity of F^l'^^*'' (i, cj) with respect to lo 



Let us define tl — j^- From (32) we have, for all k: 



ds 



tl c 



ds\ 



Without loss of generality we may consider that the external current is uni- 



formly bounded: 3i+ < +oo such that, Vj, Vs, 



< i+. Then, 



(ext) , 



< 



'El 


i+1 
1 


/* 




c 


J —oo 



< 



'El 


i+] 
V — 


f 




c 


<J —oo 



< 



+ 



El^i_ 
Tl c 



t N 

^varm\rkj{s,t, .)] + ||r(s,<, .)\\var„^[xj]{s,uj)]ds, 



from the bounds on warm[rfcj(s,i, .)xjj(s,w)] found in the previous section. 
Therefore, from (55), 



< 



Tl c 



2-Y,GM — )e-^ I (t-s)e(''-'^)(*-^)ds+ / \\T{s,t, .)\\var^[xjj{s,u;)] ds 



N 



1=1 



Using the same arguments as in the previous section we get that V^'^'^*^^ is con- 
tinuous with respect to the raster uj with an exponentially fast convergence of 
its variation. 



7.2.5 Continuity of with respect to w 

We have: 



varraQ{t, .) < 



^2 pt 



llr(s,<,w)x(s,a;)r^(s,t,a;)|| ds 



a 



2 /-t 
B 



< '-f I 2\\T{s,t,.)\\var^[T{s,t,.)]\\ds 



,2 ft 



< 



2g2(fc-i^)(t-.)(g^ 



i[G(t.)l 



(t-s) _ 



l)ds, 
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The last inequality comes from the fact that the argument of the integral is 
always positive. Then varmQ{t, .) converges to as m — > +00 from the same 
arguments giving the continuity of V^^'^^'* . 

7.3 Uniqueness of the Gibbs measure 

The uniqueness of the Gibbs measure associated with a system of transition 
probabilities P [w(n) | w"^^ ] is given by the following theorem [25j . 

Theorem 4 Let: 

m(P) = inf inf P \uj(n) \ lo1z} 1 

and 

v(P) ^ sup varn-m'[^[uj{n) | . ]]. 

If m(F) > and v(P) < 00, then there exist at most one Gibbs measure consis- 
tent with it. 

We have not been able to prove that those conditions hold in general. We 
only give here the following argument. Assume that the transition probabilities 

are differentiable ^ functions of x{n,uj) =^ ( F^'*' (n, w), Q(n, w))), and that the 
derivative is upper bounded in norm. Then, the variation of P [ ^(n) | w"^^ ] is 
bounded by the variation of x(n,uj) and this variation is exponentially decaying 
(this follows from the proof of continuity of V'^'^'' {n,uj) and Q{n,uj)). As a 
consequence, v{F) < 00. 

To ensure that m(P) > 0, one needs to slightly modify the model, replacing 
the usual reset definition where the membrane potential stays during the re- 
fractory period, by a reset to a random value, as done in |12j . 
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